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ABSTRACT 



A method to calculate the statistical properties of microlensing light curves is 
developed. The approach follows works by Deguchi & Watson, Seitz & Schneider 
and Neindorf, attempting to clarify the ideas involved and techniques used in the 
calculations. The method is then modified to include scattering by multiple lensing 
planes along the line of sight and transition to a continuous limit of this treatment 
for average quantities is performed leading to a Fokker-Planck type equation. The 
equation is solved for a particular model of the random star field and microlensing 
effect on the flux temporal variability is extracted. Applications in astrophysically 
relevant situations are discussed. 
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1 INTRODUCTION 

Gravitational microlensing of distant sources has been de- 
veloped as a tool to study dark obj ects for more than 40 
years now sin ce p ioneering work s bv iRefsdal |l964, 196q). 
iLiebeJ (|l96^ and lBvalkol (Il969l) (the idea can be traced, of 
course, to even earlier works as listed in the Introduction of 
Schneider, Ehlers & Falco (1992) monograph) and has at- 
tracted considerable interest since the effect was discovered 
in the images of the macrolensed qua sar QSO2237+0305 
ilrwin et al.lll989l : ICorrigan et al.lll99H) . In the years lead- 
ing to the discovery - and even more after it was an- 
nounced - the researchers working in the field turned their 
attention to the systematic study of the statistical proper- 
ties of microlensing as it was soon realised that the num- 
ber of microlenses significantly contributing to fiux vari- 
ations is normally large and therefore individual events 
in the light curves can r arely be unambiguously modelled 
JSchneider fc Weis jlT987ll . 

In most cases of interest, gravitational microlensing can 
be studied in the geometric optics approximation together 
with the assumption that the deflectors' gravitational field 
is weak and static. The progress in this field of research has 
been made along a number of different complimentary di- 
rections - both in terms of the questions answers to which 
were sought and the methods used in this search. On the one 
hand, many important result were obtained for the proba- 
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bility distribution function of the microlensing magnification 
value p. For example, on the analyti£_sid£^the large_£t_be- 
havio ur of this function was found (|Peacockll986t ISchneiden 
Il987l) . However, only numerical methods - most noticeably, 
the ray-shooting (Kayser, Refsdal & Stabell 1986; Schneider 
& Weiss 1987; Wambsgans s. Witt and Schneider 1 992) a nd 
image-tracking techniques JLewis et al.1.1993. ; Wittlll993r) - 
turned out to be capable of the determination of this func- 
tion over most of its domain. These focused on extracting the 
magnification PDF from many realizations of the magnifica- 
tion value as a function of the source position with respect 
to the microlenses. Starting from a simple model of static 
shearless random microlens field this approach has been fol- 
lowed to include the external shea r an d even the random 
proper motions of the microlenses IIKundic fc Wambsgansa 
ll993l:IWambsganss fc Kundidll995tlSchramm et alJll993l) . 

Perhaps less attention has been paid to a related ques- 
tion of correlations between the magnification values for a 
source placed in different positions in the microlens field. 
The autocorrelation function of this kind has been studied 
numerically using large representative samples of the light 
curves of sources tracked through magnification maps for 
point-like and extended sources (Wambsganss, Paczyhski 
& Katz 1990; Lewis & Irwin 1996). On the analytic side, 
Deguchi & Watson (1987) first calculated the variance of the 
magnification value for a Gaussian source and then Seitz and 
Schneider (1994) presented a general theory for the calcula- 
tion of the autocorrelation function of a static random field 
of stars. However, although they developed a comprehen- 
sive theory to calculate the correlated deflection probability 
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distribution used in calculating the autocorrelation function 
its practical application turned out to be rather difficult for 
technical reasons. A further major step was taken by Nein- 
dorf (2003) who found a simpler analytic expression for a 
restriction of correlated deflection probability Fourier trans- 
form sufficient to calculate the flux autocorrelation function 
for a circularly symmetric Gaussian source. 

One conspicuous distinction between the two classes of 
works mentioned above is that only the former need to solve 
the lens equation. The reason here is that the quantity they 
are concerned with - the magniflcation factor /x or the flux 
associated with a given source - is speciflc to a particular 
point in the observer's plane while the initial data to cal- 
culate it is contained in the intensity field. Physically, the 
intensity represents the density of the photons in the phase 
space spanned by their coordinates and linear momentum 
projections - hence, the use of 'phase space' in the title of 
our study. 

Viewed in this way, solving the lens equation amounts 
to 'deprojecting' the fiux distribution onto the intensity 
field; this mapping is normally one-to-many and research 
thus focuses on the deformations of the infinitesimal light 
beams and multiplicity of solutions bringing about the caus- 
tic structure. On the contrary, the flux autocorrelation func- 
tion is defined as the average of various moments over the 
observer's plane. Since every light ray eventually reaches 
some point in this plane one deals with the transformation 
of the intensity field as a whole; the conclusions are drawn 
from this transformed data. There is no need to consider 
caustics in this case since there are no caustics in the phase 
space. 

The purpose of the present paper is firstly, to spell out 
more clearly how the intensity field changes under the in- 
fiuence of gravitational microlensing and how it can then 
be used to calculate the quantities of observational inter- 
est. We try to identify explicitly which physical assumptions 
are responsible for a number of properties of this transfor- 
mation which make calculations easier and more physically 
meaningful. We show that the conservation of the surface 
brightness in gravitational lensing - which is an analogue 
of the statistical Liouville theorem - leads to the linearity 
of the transformation of the intensity n-point probability 
distribution functions and their moments. This leads to a 
slightly different description of the gravitational microlens- 
ing which is, perhaps, most useful in analytic studies, since 
the dimensionality of the problem grows rapidly with the 
transition from the configuration to phase space making it 
less amenable to numerical investigation. 

We highlight how an obvious distinction between the 
free propagation of light rays and their defiection on the 
lens plane refiects in the structure of operators performing 
the linear transformation mentioned above and in particular, 
how certain statistical symmetries of the defiector field select 
a basis in which the defiection operator assumes diagonal 
form. Dependence on the time coordinate is also brought in. 

Another advantage of the phase space description is 
that virtually no effort is needed to incorporate multiple lens 
planes along the line of sight. Perhaps less indisputable but 
yet more appealing is the transition to a continuous limit 
instead of many lens planes. We argue that in doing so one 
does not have to abandon the thin lens approximation when 
the average quantities such as the autocorrelation function 



are considered and their transformations can indeed be ap- 
proximated as continuous; this is similar to a macroscopic 
description of the scattering phenomena although their in- 
ternal mechanics is microscopic. After an argument in favour 
of the suggested approximation we obtain an equation which 
describes the intensity moments propagation in space filled 
with randomly placed individual deflectors. 

The rest of the paper is more illustrative and rather 
technical in scope. We consider the propagation of the first 
order moments which are of very limited observational in- 
terest on their own but will be the basis for the study of the 
autocorrelation function propagation in the companion pa- 
per. We calculate the Fourier transform of the defiection an- 
gle and potential time delay probability distribution density 
and try to explicitly identify the physical meaning of var- 
ious results and techniques used. The expression obtained 
turns out to be weakly dependent on the time coordinate 
thus reducing to a classic Katz et al. (1986) result. In ad- 
dition, we find that this function is dominated by the first 
nontrivial term over the entire domain of interest and illus- 
trate how the coefflcient in front of it can be obtained in a 
simple alternative way. 

Then the moments propagation equation with the ap- 
proximate potential is solved in the case when parameters of 
the deflectors are assumed constant along the line of sight. 
The basis of the solutions and the associated eigenvalues are 
found which can be used to solve the Cauchy problem for the 
equation. A particular example of an isotropic source with 
an arbitrary spatial proflle is considered and a fllter function 
describing different harmonic modes damping is calculated. 
The last section discusses the results and their connection 
to observationally accessible variability timescales. 

Throughout the paper we consider flat and static geom- 
etry and the results therefore can only describe local Uni- 
verse z < 1\ however, this can be generalized to different 
universe models although the resulting differential equations 
will no longer be autonom ous. 

Paper II of this series iJTuntsov fc Lewisll200(ll presents 
the calculation of how the flux autocorrelation function and 
power spectra are affected by gravitational lensing following 
the approach present in this study. The results are applied 
to modulation of pulsars' emission, coherence estimates for 
quasi-periodic oscillations in dwarf novae and low-mass X- 
ray binaries and the interpretation of power spectra of active 
galactic nuclei. Microlensing- induced scatter in brightness of 
distant supernovae is also considered. 



2 MODEL 

In this section, we introduce the model used to character- 
ize observational properties of lensed objects. The idea is 
to follow the propagation of the statistical properties of the 
intensity fleld and its moments from one lensing plane to 
another. Later we will attempt a transition to the contin- 
uous limit of this treatment. The geometric optics approx- 
imation is used and the angles involved are assumed small 
throughout; this is ty pical for studies of mo st aspects of the 
gravitational lensing JSchneider et al.lll992h . 

The advantage of using the intensity - that is, the 
amount of radiation energy traversing a unit area per unit 
time per unit solid angle ~ is that this quantity is directly 
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Figure 1. A schematic view of a light ray, deflected by a number 
of parallel 'lensing planes' and travelling freely between them. 
The coordinates and momenta of the ray at every plane are given 
by r and ex, as shown for the first plane, the positions which 
correspond to the same coordinate r in different planes are those 
connected by the rays perpendicular to the planes. The planes 
are distinguished by the depth parameter D, its axis is directed 
towards the observer. 



connected to the observable quantities such as flux and, 
by virtue of the Liouville's theorem, is conserved whenever 
no absorption is involved and photon frequencies are not 
changed. This includes the case of gravitational microlens- 
ing where such conservation is referred to as that of the 
surface brightness. 

The formal connection between the intensity and the 
phase space density can be made by noticing that for parax- 
ial rays the variable component of the momentum of a pho- 
ton is proportional to the direction of the hght ray, its energy 
is not changed between the emission and observation points 
(even if this is the case, the correction is simply the fre- 
quency ratio cubed), and the dependence of the line-of-sight 
coordinate on time is always one-to-one for such rays. 



2.1 Statistical description and the ensemble 

Now assume that we define a number of surfaces between the 
source and the observer in such a way that the full knowledge 
of the intensity field on one of them is sufficient to calculate 
it on the next one. 

Th e simplest choice is a numbe r of the so called lensing: 
plane s JBlandford fc Naravan||l986| : [Nottale fc ChauvinearJ 






ll986l:lKovneJl987l:IJaroszvnskilll989r) orthogonal to the line 
connecting the source and the observer. The planes can be 
labelled with parameter D starting with D = at the source 
(Figure^. At this stage, we do not include any relativity in 
the analysis and therefore our scope has obvious limitations. 

For the purpose given above, a sufficient way to char- 
acterize the intensity field on a certain surface is to specify 
the values of the radiation intensity for every point r on the 
surface along every direction a in the forward direction and, 
possibly, at every instant t, I{r,a,t). Note that in the light 
of the paragraph just before the beginning of this subsec- 
tion, these coordinates fully parameterize the phase space 
of the photon dynamics M B {r,a,t). 

Often such a complete description of the intensity field 
is not available but it is possible to use a statistical descrip- 
tion instead - e.g., the n-point probability distribution func- 



tion Pn{{I}n, {r, a, t}) defined as the probability to observe 
intensity values less than /i, ..., /„ at the respective coor- 
dinates (ri, oil, ti), ..., (r„, a„, t„). (We denote these arrays 
of I and (r, a,t) as {I}n and {r,a,t}„ to avoid unwieldily 
long notation.) 

The probabilities here can be defined as a fraction of 
appropriate realizations of the intensity field within a certain 
ensemble. One can think of two major ways to define such 
an ensemble in the context of microlensing. 

Firstly, we can consider a number of identical at D = 
realizations of the intensity field in faraway patches of the 
sky. As the light rays propagate towards the observer they 
encounter different realizations of microlenses and the in- 
tensity field realizations are subject to individual transfor- 
mations on their way. As a result, J(r, a,t; D) is a random 
function with parameter D. Our aim is to understand how 
its distribution functions P„ change with D. 

Alternatively, one can consider a single source and as- 
sume that on a sufficiently long time scale, the light rays 
emitted by the source at separate time moments encounter 
independent realizations of defiectors on their way to the 
observer. In this case, the elements of the ensemble can be 
obtained by observing the intensity field at different times. 

An important point is that the time separation between 
such observations should be long enough that the defiectors' 
configurations at these moments can be considered indepen- 
dent of each other. In practice this means that whatever 
our observable is, the time separation between individual 
measurements used to obtain the average of this observable 
must be much longer than the typical correlation time for 
this quantity fluctuation due to an individual lensing plane. 
The second alternative can only be sensibly applied to a 
statistically stationary source. 

It is normally assumed that the knowledge of n-point 
distribution functions of all orders provides the most com- 
plete statistical description of the intensity random held. 
However, just a few first ones and their moments can be 
useful in connecting statistical properties of the deflectors 
to the information that is available through observations. 

For instance, since the angular coordinate is normally 
unresolved in gravitational microlensing observations, the 
quantity of interest is the flux F - the integral of the in- 
tensity over all possible directions. Defined in this way, it 
represents a linear functional on the intensity field 



F[I]{r,t-D)= / d'aI(r,a,t;D), 



(1) 



and provided that the values of the latter are bound, it is 
easy to see that the expectation value of the flux is the same 
functional on the expectation value of the intensity: 

EF[J](r,i;D) = d^aE I{r,a,t;D) = F[EI]{r,t;D). (2) 

This is useful when we have a method to calculate the ex- 
pectation value of I{r,a,t) developed below. 

It is not difficult to devise a similar linear device to ob- 
tain the expectation value of the product of fluxes at two 
different points - i.e., the flux autocorrelation function. One 
would need, however, to deflne it on the intensity autocor- 
relation fleld to preserve linearity and thus the simple way 
to calculate the expectation value. 
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2.2 Linearity of the intensity field transformation 

Once an n-point distribution function P„{{I}„,{r,a,t}„) 
has been specified on a certain surface, its subsequent prop- 
agation is described by a linear relation. 

Indeed, assuming there is no absorption, the intensity 
is conserved along the ray, and the process we are dealing 
with is just a sort of mixing - the probability P^ to observe 
a certain configuration {I}n of the intensity field at points 
{r, a, ty„ on the next surface is just that P„ on the initial 
surface at points {r, a, i}„ connected to {r, a, t}'„ by the 
fight rays: 

P,; ({/}„, {r, a, t};) = P„ ({/}„, {r, a, t}„). (3) 

The case of interest is when the mapping {r,a,t}n -^ 
{r, a, t}'n is itself of random nature measured by the proba- 
bility density p„ ({r, a, i}„ — » {r, a, t}'„) for light rays start- 
ing at {r, a, i}„ on the first surface to get to {r, a, t}',^ on 
the next one. Then, the transformation law for P„ is given 
simply by the total probability formula: 

P; ({/}„, {r, a, iK)= (4) 

d{r,a,t}„p„ ({r,a,i}„ -^ {r,a,t}'„) P„ ({/}„, {r, a, t}„) 



where d{r,a,t}„ = Y[ d^rid^atdti was used for brevity. 
1=1 

Strictly, this expression is only valid when no point 
{r, a, t}'„ is shadowed - that is, all possible mappings in the 
backward direction lead to a certain point {r, a, t}„ at the 
preceding plane; this is normally warranted by the symme- 
try of p„ and its normalization. Even if this condition is not 
satisfied, the results are still valid for monomial moments of 
the intensity. 

The 'transition probabilities' p„ are independent of I 
and therefore the same transformation law is valid for the 
probability densities d"P/d{/}n and its various intensity 
moments - e.g., the intensity n-point correlation function: 



M„{{r,a,t}„) 



h 



/„d"P„ ({/}„, {r, a, i}„) (5) 



Therefore, the set of all such intensity moment func- 
tions on A^" naturally acquires the strutcure of a linear 
space V(A1"); this is useful because the properties of the 
transformations applied to the intensity moments often sug- 
gest particularly convenient bases in which to represent Mn ■ 

2.3 Deflection and focusing 

From the intensity field point of view, light propagation in 
the presence of gravitational lensing involves two rather dis- 
tinct kinds of transformations. 

The curved part of the light ray trajectory is normally 
very small compared to other important length scales of the 
prob lem and one adopts the so -called thin lens approxima- 
tion JPvne fc Birkinshawlll99dl - the light rays are assumed 
to be simply broken on the plane by a certain angle depend- 
ing on the deflector configuration, they also receive certain 
time delays; we will call this deflection and denote the appro- 
priate linear operator $. Between such breakages, the light 
rays propagate freely, converging or diverging depending on 
the angle between them; this will be termed focusing, F. 

In the thin lens approximation (sometimes termed just 



the gravitational lensing approximation), the spatial coor- 
dinates of the light rays do not change when a deflection 
occurs, only the angles are modified and time delays are ex- 
perienced. In focusing, one has an opposite situation. This 
refiects in the structure of the kernels (p„ in notation of I^J) 
of linear operators which describe the two transformations. 
Deflection should be proportional to identity in its spatial 
part, i.e. <I> ~ 5 ({r}„, {r}^). The same applies to the angu- 
lar part of the focusing: P ~ 5 ({a}n, {q}Jj). 

Another major difference between deflection and focus- 
ing is that only the former can be of random nature. The 
result of focusing over some distance D is always the change 
of r to r + DoL and off to t + D{l + o? /2)/c and since the val- 
ues of the intensity themselves do not change, the focusing 
operator F{D) reads 

P(i3)M({r,Q,f}„)= (6) 

M {{r - aD, a, f - (1 + a^/2)D/c}„) , 

or, for sufficiently smooth functions M{{r,a,t}), 

F (D) = exp \-D Y^ (a, • 9,. + c'^ (l + a?/2) dt,)] (7) 

with i running through all rays from 1 to n. 

On the contrary, defiection may contain both the de- 
terministic and random parts; because of the linearity, they 
can be applied separately. A simple but important property 
of the deflection comes from the assumption that the angles 
a are small. Taken with the thin lens approximation, this 
means that all rays are deflected by the same amount re- 
gardless of incidence angle a. Therefore the probability den- 
sity pn for $ can depend on the difference between a and 
a' only. This symmetry reflects in the fact that $ assumes a 
diagonal form in an eigenbasis of this transformation group 
- e.g., Fourier harmonics exp[— ir ■ a], where it is especially 
simple to see what the coefficient in front of the 5-function 
is. 

For the deterministic part - deflection /3ti(r, t) and time 
delay Sd(r, t) due to a known deflector - a' = a + fid, t' = 
f + Sd and therefore 



c&d ({r, T, f}„, {r, T, i};,) = 5({f}„ + {s}„, {t};o 



(8) 



[' 



X exp |i > Ti ■ /3d{ri,ti 



5(Mn,{rK)5(W„,{rK). 



In the case of an unknown deflector, one finds the average 
result of deflection using the expectation value of the expo- 
nential factor (E (I>M) = (eI>) M): 



*(Wn, Mn, {i}n, W;) ~ Eexp 



['E' 



l^{r^,t^] 



(9) 



- i. e., f3 characteristic function, Fourier transform of its 
distribution function density p({/3}; {r}, {t}). 

Note that this trick does not work for time delays: since 
the addend s to f depends on t itself, exp[iw s]-like factor 
would not come out of the integral as a factor when mak- 
ing the Fourier transform in t. However, when the expecta- 
tion value Eexp[iajs] is independent of ti for all i — l,n, 
<1> assumes a similar form in its temporal part. This is the 
case for a statistically stationary deflector when only one-ray 
statistics are considered. This problem, n = 1, is considered 
below. 

Of our prime interest are deflections due to a priori 
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unknown distribution of (niicro)lenses and we will ignore 
the deterministic part for the rest of this study. 



2.4 Multiple lens planes and the continuous limit 

The above approach, in slightly different expression and with 
a different emphasis, has been used to study the deflection 
statistics by a single deflecting plane (e.g., Katz ot al. 1986; 
lNeindorl2003^ . Little change is needed to apply it to a num- 
ber of different planes when the light rays are significantly 
deflected a number of times on their way from the source to 
the observer - one applies the deflection and focusing trans- 
formations a number of times with appropriate parameters 
thus implementing the multiple lens plane approximation 
within this approach: 

Mo ^ F[Dok)^kF{Dk,k~i)-F{D2i)^iF[D^s)Ms, (10) 

where Mo and Ms are the average intensity at the observer 
and the source and Djj-i is the distance between the j'-th 
and {j — l)-th lensing planes. 

Focusing operator is diagonal in spatial wavevector k, 
incidence angle a and temporal frequency uj representation, 
as can be seen from Q. If one disregards potential time 
delays, the deflection operator assumes diagonal form in r, r, 
oj basis. The conversion between the two representations can 
be efficiently performed numerically which suggests using 
the formulation presented for numerical studies of multiple 
lens planes case. 

Our focus, however, is different. When one deals with a 
3-D distribution of lenses such as a cosmological population 
of primordial black holes, CDM haloes or dense gas clouds, 
this somewhat arbitrary 'slicing' of 3-D space into a number 
of lensing planes seems a little unnatural. In addition, since 
only limited information about the distribution of such de- 
flectors is normally available (or assumed) , the main interest 
is general behaviour of various observables and their aver- 
ages as the light propagates through the Universe. We are 
therefore naturally lead to a three-dimensional problem. 

The complexity of the full-scale problem is enor- 
mous and only few attempts to attack it directly have 
been made. There are two main reasons for this. First 
is the need to abandon the thin lens approximation 
- one cannot describe a gra dual change in light ray 
momenta within this model JPyne fc Birkinshawl 1199' 



lYoshida. Nakamura fc Omotell200a) . The second is also re- 
lated to the thin lens approximation - when the deflection 
probabilities are computed, the correlations between the de- 
flectors belonging to different planes would have to be con- 
sidered. 

However, when it comes to the statistical average of the 
quantities of interest, the thin lens approximation can still be 
valid for two reasons. Firstly, since the incidence angles are 
very small, deflections due to different individual lenses sim- 
ply add up. Secondly, in astrophysical applications, a typical 
scale over which moments (|SJ| change appreciably is much 
larger than the light ray curvature radii. Therefore the sta- 
tistical averages can be taken on 'statistically' infinitesimal 
dphya^ scales which are much larger than the scale relevant 
to elementary deflections. 

Since the deflection angles are additive, averaging of 
dphys?E> reduces to the use of the thin lens approximation. 



while the commutator [dphya$, -F(dphya-C))] , which quanti- 
fies the error introduced, is of the second order in dphys^- 
From now on, when some changes in the intensity distribu- 
tion or moments are mentioned, the averages over the same 
physically small distances will be actually meant; the sub- 
script 'phys' in differential will also be dropped. 

The above observation can also be used when com- 
puting the deflection probabilities - whenever the scale of 
correlations between individual lenses is small compared to 
dD, their distribution can be considered effectively two- 
dimensional. The model distribution used in the next sec- 
tion assumes that no correlation at all is present and all 
individual deflectors are distributed independently. We are 
not aware of any similar calculations which would take the 
correlations into account. 

We have thus put some additional meaning in the thin 
lens approximation and can now obtain a differential equa- 
tion to describe the continuous propagation of the intensity 
field moments. Indeed, the changes in the intensity field av- 
eraged over the ensemble are small on a small run dD. There- 
fore we can expand F and "1> operators in series of dD. The 
zeroth order term in this series is identity, and therefore, to 
the first order in dD, 



dM ,dF d-l- , , ^ 



(11) 



- a result of a typical Fokker-Planck structure dtp = Lp. 

For independently distributed deflectors, the only ex- 
tensive parameter describing the distribution is the total 
number of the deflectors TV. Since the two otherwise similar 
transformations with parameters A^i and N2 applied succes- 
sively should produce the same result as a single transforma- 
tion with the pooled parameter: l>(Ai'i)<l>(A''2) = <b{Ni+N2), 
there exists an operator if> such that: 



$ {N) = eyip[Nip\ 



(12) 



and since the number of deflectors in a layer of depth dD is 
proportional to dD, one has 



d$ 
dD 



^ = i^m 



dexp[A''i^] 



diV 



i/|l£||(^ 



(13) 



where u is the space number density of the deflectors, ||£|| 
is the area of the lensing plane C and the operator if> only 
depends on the individual properties of deflectors. 

In the case of the first order moments (or the 1-point 
intensity distribution function) the equation simplifies sig- 
nificantly when one assumes (statistical) stationarity of the 
deflectors. In this case, the expectation Eexp[i(T ■ (3 -\-u a)\ 
is independent of t and ip assumes diagonal form in r, r, lo 
basis. 

There is a certain trade-off between the simplest forms 
of the focusing and deflection - one could instead leave the 
dependence on angles to make the focusing algebraic but 
that would leave the deflection an integral operator which 
is harder to deal with. In addition, many properties of the 
solution will be dictated by the form of the focusing which 
does not depend on the deflector properties and distribution; 
we therefore opt for r,T,uj coordinates. 
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Changing to this basis and evaluating AF/AD at D = 
using ||7|l one arrives at a Schrodinger-Uke equation^; 

'idnM = ^ dlM - dr ■ d^M + WM (14) 

with a complex potential iV{v,T,uj) = id"I>(r, T,Ci;)/dD: 

V{r,T,iu) = u\\C\\^ln{ENe^p[i{T-^{r)+u;s{r))]}{15) 

where the subscript N at the expectation symbol E means 
that functions /3 and s are computed for the distribution 
with the (average) number of deflectors A'^; the derivative is 
evaluated at N = 0. 



3 DEFLECTION PROBABILITIES FOR ONE 
RAY 

We e mploy the Markov-Ho ltzmark-Chandrasekhar method 
(e.g., IChandrasekharl Il943r to calculate the deflection po- 
tential V due to randomly and independently distributed 
deflectors. The essence of this approach is in the fact that 
the deflection operator for one ray is diagonal in r, r, oj ba- 
sis for a statistically stationary deflecting plane, as explained 
above. 

Therefore, transformation due to exactly A'^ deflectors is 
just the N-th power of transformation due to one deflector: 



p{N;r,T,uj) = (p(l;r,T,tj)) 



(16) 



This result has been used by a number o f authors in gravita- 
tional lensing studie s (Katz et al. 1986; IPeeuchi fc WatsonI 
ll987l:lNeindoril2003l) . 

Since there is no control over the number of deflectors A'^ 
in the layer, and especially since we consider a 'physically' 
infinitesimally thin layer dD, the number TV is allowed to 
fluctuate. Assuming Poissonian distribution for thi s quantity 
with average N, the average transformation reads (jNeindorj 
l2003l): 



Pjv (r, T, uj) = exp[7V (p (1; r, r, uj) - 1)] . (17) 

and using 11511 one has for the deflection potential; 

V{r,T,iu) = u\\C\\{p{l;r,r,iu) ~ 1) (18) 

= i/| iril {E exp [i (r ■ /3, + (X) sO] - 1} > 

where f3i and Si are the deflection angle and the potential 
time delay due to a single deflector which depend on its 
individual properties only; the expectation is taken as the 
average with respect to all these properties. 

An individual deflector can be characterized by the 
above two functions /3i (r — r') and Si (r — r') that, in the 
deterministic case, describe the deflection angle and poten- 
tial time delay observed at r due to a deflector placed at 
r'^. They depend on other deflector parameters as well, 
for instance, they must be proportional to its mass m. In 
the simplest case of a point-like mass f3{p) = —mp/p^ and 
s{p) — —m/c\np/po for p greater than a few gravitational 
radii m/2, po is a constant; it is convenient to choose it in 



^ We drop u>/c term which describes the fiducial light-travel time 
and is ehminated by the transformation M — > M exp[iijjD/c]. 
^ In gravitational lensing f3i is, up to a factor, a gradient of Si, 
but this will not be used here. 



such a way that the average potential time delay is zero 
(Schneider et al. 1992). 

In this study we restrict ourselves to the simplest case of 
a point-like mass; the only parameter m is flxed so that the 
expectation value is calculated as the average with respect 
to the deflector position in the lensing plane £,: 



E/: 



(r')/(r') 



(19) 



r'e£ 



where is the distribution function density of r'. 

It is most natural to assume that individual deflectors 
are uniformly distributed in C: 4>{p) = ||/I||~'^, so that the 
resulting potential is independent of r in the limit of infinite 
lensing plane C. 

However, the radial integration in 11911 for r~^ exponent 
diverges both at infinity and at zero. The latter limit poses 
no problem because the integration should be truncated at 
some finite distance from the deflector e = 0(771) - the light 
rays cannot pass through if shot too close to the lens. As 
for the outer integration limit, no satisfactory explanation 
appears in the literature for putting a certain cut-off. The 
discussion normally revolves around the fact that all real 
distributions have some physical scale associated with them 
- e.g., the size of the galaxy the deflectors belong to. 

However, for a flnite lens plane, the boundary effects 
start to matter, and what is difficult to secure is p's in- 
dependence of r. Standard derivations (Katz et al. 1986; 
Schneider et al. 1992) rely heavily on observation point r 
being the symmetry centre of the distribution and it is easy 
to check that once this symmetry is broken p starts to de- 
pend on r; moreover, this dependence does not vanish in the 
limit of an infinite jC. 

This problem appears to be similar to standard cosmo- 
logical paradoxes involving the infinite distribution of mass 
and light in the Universe and the absence of the satisfactory 
solution to it might also be an indication that the consistent 
solution can only be found within fully relativistic approach. 
However, as the latter is not available at the moment, we 
will use a standard way out: to imagine that £ is symmetric 
enough to make r dependence vanish locally, calculate p for 
such a distribution and then choose the lens spatial dimen- 
sion in a way appropriate for the problem. This is expected 
to capture the general behaviour of the intensity moments. 
However, one should keep in mind limitations inherent in 
this derivation. 

Armed with this knowledge, we calculate p{t,uj) = 
p{l;r,T,uj) for a defiector distributed in a disc of radius 
R for observation point r = 0. That is, we choose the dis- 
tribution 






e^ ps; 7? 



(20) 



and zero otherwise; e <C -R is assumed, so that <j}{p) is 
still normalized. We use I3i{p) = —mp/p^ and Si(p) = 
—m/c\np/ pQ, Po parameter will be chosen shortly. 

To calculate p{t, uj) we expand the exponent in its Tay- 
lor series and then integrate it term by term: 

p(t, u) = E{ 1-l-i {t-I3,+u; sO- 2 [r ■ |3^ + uj s,f + ...}. (21) 

The resulting series converges for (p defined by 12UL and 
therefore equals the initial integral. 
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The actual integration is done in Appendix^ and the 
result is presented below. However, the expansion above is 
quite useful in clarifying the physical meaning of the po- 
tential function V , proportional to p — 1. The first term is 
the normalization, integral of p over all /3 and therefore V 
expansion starts from zero. 

The second term contains average r • (3 and u) a. The 
former has to be zero because of the symmetry: all the di- 
rections of /3 are equally likely; for the same reason all other 
terms containing odd powers of r vanish. As for the average 
potential time delay Es, we can use our freedom in choosing 
po in Si{p) to make the average potential time delay zero as 
there is no way we can measure this delay directly. 

Therefore, the expansion of V starts with terms pro- 
portional to the squares of r and uj, the cross-term in the 
square of the sum vanishes because it contains an odd power 
of /3. In fact, the quadratic term in r, proportional to the 
dispersion of /3 dominates V , the correction is a slowly vary- 
ing function of order ln~^ A, where A = R/m, and one can 
write, for lu = Q: 



P[T) 



^E/3^ 



(22) 



where an extra 1/2 is the average value of cos^(/3,r). An 
analogous interpretation for uj ^ can hardly be provided 
because the deflection angle and potential time delay are 
strongly correlated - in fact, /3 uniquely determines s. 

The average square of /3 (which is also its dispersion) 
can be calculated in a simple way by inverting /3(p) to obtain 
the distribution density of the deflection angle 



p(/3) = 0(p(/3)) 



D{p) 



Dm 



\C\ 



m 



Observing that /3 ranges from /3min ~ m/R to /3„ 
uses the above distribution to calculate 



Efi^ 



2tt- 



\C\ 



■In A 



and at UJ ^ one has 



Vir) 



■In A. 



(23) 



1 one 



(24) 



(25) 



More accurate calculation, performed in Appendix IXI 
results in the following deflection potential: 



V{t,Uj) — TTV 



7?^A/.(-!^,A)+mV„(-!^,A)](26) 



'^ 2 2,. 

— vm T InA 
2 



1 exp(ia;/t^o) ^ 



exp{iuj/ujo) 

ibj/ujo ' InA 

and A/cj are deflned in Appendix 1X1 and 



Functions /r, fu 
Wo = c/(m In A). 

The flrst line in the above expression, which is inde- 
pendent of T, describes some additional spread in light ray 
arrival times due to the variance of potential time delays 
within the ensemble. It is easily taken into account and will 
be discussed in greater detail when solving the equation l|14|l 
in the next section. 

For the second line, which represents the r-dependent 
part V' of the deflection potential V, one can notice that 
for a solar-mass lens and In A ~ 10 — 10^ , the characteristic 
frequency wo ~ 5 x (10 — 10 ) s~^. This is much higher than 
typical observationally measurable frequencies in astronomy, 
and in the leading order, one obtains 



V'{t) 



^ 2 2,. 

— vm r InA 
2 



1 



InA 



U 



(27) 



This approximation for V' is independent of the frequency 
Lu. Physically, neglecting this dependence means ignoring the 
potential time delays altogether, as can be seen from 11511 . 

This can indeed be done due to the logarithmic depen- 
dence of the potential time delay s on the impact parameter 
p: p ^ exp[— s/2tm], tm ~ m/2c is the gravitational radius 
light crossing time. The average time delay is unobservable 
and we remove it by choosing appropriate po; at the same 
time, close encounters are rare for realistic deflector den- 
sities, and one can safely neglect all deflectors except the 
dominant, closest one. The distance to this deflector p in 
the layer of depth d is distributed according to the Poisson's 
law: P{p) = 1 — exp[—Tvdiyp^] and therefore potential time 
delay distribution density is 



dPjpjs)) ^^ nudpl 

ds tm 



exp 



1 2 -s/tm U 

'TTVdpQe ' — S/tri 



(28) 



- i.e., it is exponentially suppressed for s greater than a few 

In addition, because function /t is very weakly depen- 
dent on its argument (see Appendix ^ and supressed by 
In^^ A factor in 12611 . this dependence will also be neglected 
for the time being. As a result, V' reduces to 1251 1 and the 
equation 1141 1 is formally similar to the Schrodinger equation 
for harmonic oscillator with a complex constant. 



4 SOLUTION OF THE EQUATION 

We are now ready to formulate the problem properly. We 
seek a solution of the equation 1141 with a potential given 
by 12611 and initial conditions 



M(r,T,uj)\D^o =A/o(r, T, 



(29) 



When the deflectors' distribution is uniform, the potential 
is independent of r and things simplify further by changing 
into Fourier domain in r space: 

dDM{k,T,uj) =k-c»^A/(k,T,w)-^9^Af(k,T,a;) (30) 



2c 



+ V{T,Uj)M{k,T,Lj) 



Therefore, equations for different k,(.<j-modes ^ku(''") decou- 
ple for such potential, leaving k and u) as parameters of the 
equation. This starts the procedure of variable separation 
performed in Appendix [HI 

Before this, however, let us observe one further simpli- 
fication of the problems which arises from a particular form 
of the potential. Indeed, the term proportional to the flrst 
line of 12611 . that is independent of r, can be eliminated from 
the equation by making a change of the dependent variable: 



with 



w (uj, D) = exp 



■n / AD'v {R^Af^ + rnf^) 



(31) 



(32) 



which leaves only the r-dependent part of the potential V' 
in the equation for f'. 
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The physical meaning of the transformation 1311 can be 
understood using the Fourier convolution theorem. Indeed, 
multiplication of Fourier transforms corresponds to convolv- 
ing the i-represented moment M with a 'spread function' 
w{t — t') whose Fourier transform is given by 1321 . 

Therefore, when one considers a 'synchronous' source, 
in which temporal variations only modulate a given intensity 
distribution: 



Mo{KT,i^)^A{uj)I{k,T) 



(33) 



modification 13111 applies to the amplitude factor A only. As 
a result, evolution of M is split into a non-trivial part caused 
by defiections and described by equation on ^' with potential 
V' , and the overall spread in light-ray arrival times, which 
comes from the variance in potential time delays within the 
deflector ensemble. The latter part is described by 13211 . 
which is discussed in greater detail in Appendix IHl 
The resulting problem for (^'^^ 



do^'-k-d^^'~'^^dU' + V'^' 



(34) 



is solved in Appendix^using separation of variables for the 
-D-independent distribution of microlenses under assump- 
tion that fr term can be neglected in V' , as explained in 
the previous section. 

A complete eigensystem of solutions M and correspond- 
ing _D-eigenvalues A of 13UII is built which provides the solu- 
tion to the problem presented in the beginning of this sec- 
tion for any initial conditions 1291 as a sum (integrals with 
respect to k, tj in a non-compact case): 



M{r,T,t;D) = 

^«;(a;,D)(ML,. 

knZu; 

with the basis functions 

Afk„,„(r,T,i) = 



(35) 



,Mo)^ 



'M. 



(r,T,i) 



(36) 



(27r) 



■ exp 



■ r + Lot + nx'\ k ■ 



Cc^ni(r), 



and associated eigenvalues 



Xku^ni ^-—k^+i-x/iK{cj){l + \n\+2l), (37) 

where the 'oscillator constant' k is defined in analogy with 
QM harmonic oscillator: 



k{ijj) 



TTU In A. 



(38) 



so that at low frequencies tj <C coo one has 2cV' /co ~ kt^. 

Function C,niui{T) is the eigenfunction corresponding to 
l-th radial mode of 1341 . it is defined in Appendix IHl index 
n denotes angular modes. It is convenient here to present 
the solution as a function of r (rather than a) because the 
fiux corresponds to value of intensity moment at t = 0. 

The scalar product (•, ■) used above to project the ini- 
tial conditions onto the eigenfunctions is defined in the stan- 
dard way (first factor is complex-conjugated in the integral) 
in V{M). Because the potential V' is complex- valued, Her- 
mitean conjugates Qj^^ and i!^niuj are complex conjugates of 
each other; in this way the orthogonality relation IB121 
holds, as explained in the Appendix ^ The conjugates 
^kniLj *o Mi^rii^ differ only in the use of ("^^^ instead of 



Given the rule adopted in calculating the square root of 
in, the real component of eigenvalues is always non-positive 
and is closest to zero in the most symmetric case, n = and 
/ = 0. In the approximation we use, the absolute value of 
the decrement 



ReAt^n! 



i^^ujlnA{l + \n\+2l), 



(39) 



and therefore, for a fixed n, I mode it is proportional to the 
square root of the frequency. In particular, the decrement 
is zero at zero frequency - the fact that simply reflects the 
conservation of energy, which is only redistributed between 
different light beams. 

This dependence allows one to understand the phys- 
ical meaning of this damping. Indeed, it comes from two 
approximations: first, we have neglected the potential time 
delays, which makes k oc u!~^; second, the incidence angles 
are small, and therefore geometric time delays are propor- 
tional to squares of these angles, resulting in second-order 
differential operator in t. Therefore, we are dealing with 
a partial loss of coherency in a given cj, n, I mode due to 
random geometric time delays uncorrelated between lensing 
planes. Relative importance of such delays is different for 
different n,l modes, as described by 1391 . 



5 AN ILLUSTRATION: THE CASE OF AN 
ISOTROPIC SOURCE 

Let us now consider a test case where emission from a source 
located at D = —Ds enters the deflectors-filled region of 
space D > 0. It is convenient to rewrite the initial intensity 
moment as 

„2 



Mo (r, a,i) = / r- aDs,a,t 



2c 



-Ds 



(40) 



because for an isotropic source, the function / is independent 
of its second argument, and in the paraxial approximation 
used, isotropy can be assumed to describe most of the as- 
trophysical sources. In this case, / represents just a spatial 
structure of the source and its variability /(x,s), with x, s 
being the spatial coordinate and time at the source position 
—Ds- Making Fourier transform of 14U1 yields^: 



Mo(k, 



= /(k. 



LOD, 



exp 



2ljD 



'"^ {T + kDsf 



(41) 



where /(k,(j;) is the source Fourier transform. 

Projecting this initial conditions onto the basis vectors 
and summing the resulting series at r = 0, one obtains for 
the flux: 

ic/(k,cj) w{Lj,D)exp[--^k'^{Ds+D)] 



Fk^iD) = 



■>Ds cos iiS ^i^ _|_ 



'-^Vb^' 



(42) 



This is to be compared to the propagation of flux in the 
absence of deflectors, k = 0: 



Fi.^{D)/Fk^{D)\^=o = w{uj,D)K{uj,D), 



(43) 



which defines the 'filter function' K{uj,D). It can be ex- 
pressed as a function of two parameters, D/Ds and x = 

\/TkluD/c: 

^ Regularization might be done by adding a small negative con- 
stant to iLjDs/2c, then letting it tend to zero 
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Figure 2. Minimum of tlie mode decay amplitude l^olmin ^^ ^ 
function of xq. 



K{u),D) 



sin X cos X — sin x/x 



X 1 + D/Ds 

and tends to an even simpler form 

Ko{x) 



(44) 



(45) 



smx 
as D > Ds . 

Using IIB9|l we can separate the dependence of x on the 
frequency 



V-TTi/D^m In A\/e'"/"t) - 1 



(46) 



with u!q^ = 77ilnA/c, so that the first factor xq — 
Vyri/D^mln A determines the amphtude of the mode decay 
due to lensing while the second sets the frequency scale: the 
modulus of Kq reaches its minimum 

2x0 



l^oLi„(xo) = 



ycosh 2\/2a;o — 1 
aX uj = TTUJo, and at low freqencies u <S^ ljq one has 

, 4 



^0=1- 



ISxq + Xq ( uj ^'^ 



o 



\UJO/ 



(47) 



(48) 



180 \loq/ \ujo, 

The graph of |^o ];-„;„ as a function of .tq is plotted in figure|5| 

The effect of this filter is, perhaps, easiest to illustrate in 
the case of a single (5-like in t impulse. According to the con- 
volution theorem, this will be transformed into the inverse 
Fourier transform of Ko{xo,ui). Since the filter function is 
strictly periodic in frequency, Ko{xo,t) will only be non- 
zero at t = nto, to — 2TrL0Q^ = 27r In A m/c (n G A/"); this 
discreteness, however, is clearly an artefact of the approxi- 
mation used for the deflection potential. Figure |3 presents 
the magnitudes of these peaks K„{xo) for a number of dif- 
ferent values of xq. 

One should bear in mind that this broadening is calcu- 
lated as the average over an ensemble. In relating this to a 
single realization it can be argued that since microlensing is 
known to produce many microimages of the source - each 
with its own magnification and time delay - the result in 
figure 13 should resemble an average shape of a microlensed 
5-like impulse. However, such an interpretation comes from 
outside the model and cannot be justified within it (see 




Figure 3. The effect of applying the filter Kq {xo,ui) to a single 
impulse at t = 0. The dots show, for a number of values of xq, 
the magnitudes of the resulting peaks K„; these are joined by the 
lines of constant xq for convenience purpose only as K{xo,t) = 



at t ^ nto. The time scale unit is to = Sttloi., 



: 27r In A m/c. 



IWilliams fc Wiieral997l . however, where the results of a sim- 
ilar 2-D numerical calculation look surprisingly similar in 
terms of the time scale and the shape of the average lensed 
pulse in general). A legitimate statement would be that the 
average light curve produced by a number of S{t — is)-like 
impulses all occurring at the same moment ts will be repre- 
sented by Kn (or a similar continuous function) . It is hard to 
think of a physical reason for such synchronization, though, 
and marginalizing over ts averages the effect to zero. 

Similar consideration apply to damping of any harmonic 
modes which can be present in the initial signal. The aver- 
age of any mode amplitude over the ensemble - which is the 
only quantity the approach of this paper allows to calcu- 
late consistently - is always zero unless the initial phases of 
the oscillations in all individual realizations comprising the 
ensemble are somehow synchronized. This is why we are in- 
clined to conclude that this damping is even less likely to be 
observable in the flrst order than the broadening discussed 
above. Second order moments such as the power spectra, 
which do not suffer from the absence of phase coherency are 
considered in the companion paper. 



6 OBSERVATIONAL PROSPECTS 

To determine where the effect described in the previous sec- 
tion can be important - regardless of the interpretation - let 
us consider its amplitude parameter xq. For observationally 
noticeable mode damping one need xq > 1 as seen in the 
figure 121 There are two conceivable ways to achieve this: by 
increasing v or D. The former possibility corresponds to ob- 
serving a source through an extended overdensity along the 
line of sight; restricting ourselves to a single such overden- 
sity, we would interpret D as the spatial extent of this struc- 
ture rather than the total path length towards the source. If 
the total number of the deflectors in this overdensity is flxed, 
and it is largely spherically symmetric, one has v ~ ND~^ 
and therefore 
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Xo ' 



Nm 



In A. 



(49) 



As Nm is (twice) the gravitational radius of the entire over- 
density, it is unlikely that xq reaches a value of order unity. 
For instance, it has been argued, that the SgrA* source in 
the centre of the Galaxy could be a very dense cluster of 
stellar-mass (10 Mq) black holes instead of a single super- 
massive one, although this idea is almost ce rtainlv ruled out 
jMelia fc Falckel200ll : lMeha fc Cokeill lQOgY In this - highly 
contrived - case Nm/D ~ 0.1, Xq becomes a quantity of or- 
der unity (In A = InD/m ^ 10 — 100) and gravitational 
lensing could be an important effect when interpreting vari- 
ability of objects seen behind the cluster at frequencies close 
to ujo ~ 10'^ Hz. For more realistic conditions, however, xq 
would be much lower and the effect can be ignored. 

The other way to increase xq is to probe larger, cosmo- 
logical distances along an on average uniform distribution 
of compact objects. One can rewrite xo in more relevant 
cosmological terms then: 



Xo = \ -ildmA, 



(50) 



where 0^ is the deflectors' mass density in terms of the 
critical density and H is the Hubble constant. This expres- 
sion cannot be precise since the analysis we presented in 
the previous sections was developed with stationary three- 
dimensional geometry in mind, but it should give a fair es- 
timate for moderate redshifts z <1. 

Therefore, for any appreciable cosmological population 
of compact objects 2:0 reaches a value of order unity for 
moderate z and therefore such a population would mani- 
fest itself in suppressing high-frequency variability of distant 
sources. The characteristic frequency, however, remains very 
high for reasonable values of the compact objects' masses: 
as In A ~ 100 and c/m = 1O^(M/M0)"^ Hz, leading to 
Wo ~ 1O^(M/M0) Hz. At moderate redshifts, this rapid vari- 
ability is likely beyond the reach of the present-day obser- 
vational facilities. 

The only kind of sources which are bright enough to 
be probed at this frequency level seem to be the Gamma- 
Ray Bursts, but they are not expected to contain much har- 
monic signal in their intrinsic light curves at any frequencies. 
Being the results of multiple violent collisions between the 
shock waves in the expanding shells (e.g. (Zhang & Meszaros 
|2004) . GRB light curves do not have the phase coherency 
necessary in the first-order effect considered in this paper. 
The question of how their variability is affected by lensing 
can therefore be answered by the appropriate account of the 
second order statistics - e.g., auto-correlation functions or 
power spectra, which is the subject of the companion paper. 



7 CONCLUSIONS AND DISCUSSION 

In the present paper, we have attempted to clarify the ide- 
ology used when studying how statistical properties of light 
from distant objects are affected by microlensing. The im- 
portant steps of the associated calculations which are of- 
ten omitted and the justification for various approximations 
which are normally assumed tacitly have been listed. We 
hope this will serve better understanding and attract fur- 



ther interest to the field which, over the past decades, has 
accumulated enough data to make the statistical analysis of 
light curves routine. 

In particular, we demonstrate how the conservation of 
the surface brightness leads to a particularly simple - linear 
- type of transformations when the statistical properties are 
defined in the phase space and one talks about intensity field 
rather than fiuxes. Connection of the phase space description 
with the observable fiux, however, can only come through 
the average values of the latter, and one has to increase the 
dimensionality of the problem in order to preserve linearity 
when studying higher moments of the flux distribution. This 
is why we believe that the proposed approach will be more 
useful in analytic studies rather than numerical ones. 

We distinguish between the two types of transforma- 
tions corresponding to free propagation of light rays and 
their deflection on lens planes because these have rather dif- 
ferent properties and it is of advantage to represent them in 
different bases. The former, focusing, is not stochastic and 
in a sense, serves as a means of referencing the undeflected 
light rays. Deflection, which, under the thin lens and small 
angle approximations, assumes a diagonal form in angular 
Fourier domain, can have both the stochastic and determin- 
istic parts and due to the linearity of the transformation 
these can be applied separately. 

For the stochastic part, we calculate the Fourier trans- 
form of the deflection angle and potential time delay proba- 
bility density p{t, uj) thus extending the work by Katz et al. 
(1986). The potential time delays turned to be unimportant 
in observationally interesting frequency range, however, and 
p itself can be well represented by the first non-trivial term 
of its Taylor series expansion. 

The method is also particularly suitable to study the 
case of many lens planes along the line of sight. To illustrate 
this, we take the multiple lens plane treatment to the ex- 
treme deriving an equation for the continuous propagation 
of intensity moments under the assumption that the thin 
lens approximation still applies to the average quantities. 
The equation is solved for the first order moments affected 
by randomly distributed compact objects and the results are 
presented with an emphasis on the propagation of the flux 
temporal variability. We have devoted a fair fraction of the 
paper to the solution of the equation because it will form 
the basis to study the second order moments. 

The conclusion we reach is that when astrophysically 
reasonable parameters of microlens populations are as- 
sumed, the effect on the average amplitudes of harmonic 
oscillations starts to be important at frequencies which can- 
not be effectively probed in observations at the present time, 
and the variability expected from astrophysical sources does 
not extend to those high frequencies. This is not the case for 
the second moments as we show in the second paper of this 
series for the autocorrelation function. 

The scope for further development of the method in- 
cludes, first of all, generalization to non-static and, perhaps, 
non-flat geometry which would make possible its applica- 
tion to cosmologically significant distances beyond 2: ~ 1; 
this is also likely to eliminate some awkward inconsisten- 
cies in the calculation of the deflection probability Fourier 
transform such as the use of rather arbitrary cut-off radius. 
Another important direction is to obtain the solutions for 
more accurate approximations to the exact transformation 
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potential V as well as derive an expression for this function 
for deflectors other than the point-like one. Account of cor- 
relations between deflectors position is also of great interest 
since massive objects in the real Universe are all clustered; 
deflector proper motions also need to be considered and we 
do so in the companion paper for a Gaussian velocity distri- 
bution. 

The continuous approximation presented in this paper 
is not, however, a truly three-dimensional approach as the 
thin lens approximation is still in use and it is not clear 
how 3-D clustering of deflectors can be incorporated in this 
treatment. What the true three-dimensional statistical lens- 
ing will look like, therefore, remains to be seen. 
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APPENDIX A: CALCULATION OF THE 
DEFLECTION PROBABILITY FOURIER 
TRANSFORM 

To calculate the deflection probability Fourier transform we 
expand H18|l in Taylor series according to l|21|l and then cal- 
culate the integral 1191 1 in r', using 12011 : 






(Al) 



Y^ 



de dpp{l + } - (im (^so(p) + r • f3o{p))T 



With a change x = p/R, e^ = t/R, xg = po/R, a, = —mui/c 
and b = rm/ R = egrm/ R one has 



2tt 1 

AG 



N , / de' / , V^ / , a; ,cos6'\",. , 

p(T,Lo) -1= — dxxy -7 a In b A2 

/ TT / ^ — ' n! V xo X J 



. 2mijj / , , X 

-1 / da; a; In — 

c / a^o 



27r 1 

de 



+ I — / da;a; > —; [aln b 

n! V Xo X J 



n=2 



We choose xq such as to eliminate the second line. This 
quantity represents an overall potential time delay due to a 
smooth matter distribution associated with <^(p) and given 
a certain freedom to define the time axis, we set the factor 
at muj/c to zero. This leads to 



1 



'1/2 



lo ~ exp[-- - e^lne^] Ri e ^'^ ~ 1 



(A3) 



or po ~ R/^/e ~ R. 

One then further expands the n-th power of the expres- 
sion in parentheses of 1A2II using Newton's binomial formula 



(aln — 
V a;o 



, cos6i\" 



a;o X 



fc=0 
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and p — 1 becomes 

p{t,u)~1= (A5) 

n = 2 fc=0 g 

Only fc-even terms survive integration with respect to the 
angular coordinate producing 



=fc = i / de cos'= e = 2^^^il^, k = 2i, 



(A6) 



while for all fc-odd terms H2J-1 — 0. 

In radial integration, four rather different types of inte- 
grals arise: 

1 
da; 



k = n. 



r 



. , , ^ rk _ I <ix , n-k X 

4 < A; < n, In = I ,. , m 



X Xo 



fc = 0, 



/„ = /da; In" — , 

Xo 



(A7) 



(A8) 



(A9) 



(AlO) 



Rearranging summation order in 1A5II one obtains: 

(A16) 



00 [n/2] ,^ 



2j 



n=2 j=0 



= 2 



i:^'"n©E 



{iar ,2 



(rz-2)! 



;n„n-2j 



00 r, ■ 00 

+ Z^ (2i)! ^. (n-2j)! " 

j=2 n = 2j 



2i 



Plugging expressions for I^ above into these sums produces: 



p{t,uj) - 1 = e^/^(a,e^) + A/^(a,e^) 



(A17) 



-ia In xq/e 



pialna::o/ei _ ]^ 



+ A(-.,2 



4d 



where 



/o(2yf)-i-l 



/^ W - 2 E ^[(m + l)!F = y dt^^ii^^i^L^-^ , (A18) 



,, _!,_ l+ialna;o/6^-aV21n^a;o/e^-e-'°'""°/''" 



-^In^^ 



and 

•*- — ^ n 



(A19) 



(A20) 



which will be estimated separately. Two are elementary: 

2 — n -1 2 — n 

C = " ~ « ^^-r, n ^ 4 (All) 



n-2 n-2 



and 



1 A n-l 1 1 n-1 ^x 



a;o 



Xo 



'ly 



. In" ^a;o/£^ 
n-l 



n>2. 



(A12) 



For 4 ^ fc < n the integral is strongly dominated by the 
lower integration limit as e^, ^ 1, therefore 

,„_fcln""''a;o/e2, 



/„ ^ (-1)" 



4 ^ fc < 71, (A13) 



(fe- 2)6^2 

and one can see that at this level of approximation lAlU 
is just a special case of IIA13L so the last inequality can be 
changed to an non-strict one. 

For k = there are two distinct regimes; for large n > 

n = 21na;o/e2, 



jO ^ ("l)"£z , n+i £0 



(A14) 



n + 1 ex 

while for smaller n /„ should be computed by direct inte- 
gration or using a different approximation 



^" ~ 2" + i " 



(A15) 



We will write 1° = 1° + AI° with 1° defined by the right- 
hand side of (IA14I I and A7^ becoming unimportant for n > 
n. 



should be well represented by its first ~ n{ex) terms. To the 
first order of approximation it is just — a'^a;o/8. 

Returning to original variables, one obtains (A — R/m): 



, , , m^ , / muj . \ . r ( rnuj 



A 



2 2 

-e " 



2i?2 



fr\- 



1-e-'^ 



(A21) 



^) 



from which the expression 1261 follows immediately. 

Finally, we note two approximations. First, us- 
ing <A19IA20l l one obtains for the first line of 12611 at low 
frequencies: 

V{t, uj) = -Tviy^^^ (r^ + -m'' In^* a) + V'{t, uj). (A22) 
8c-' V 3 / 

Second, for large negative arguments fr has logarithmic 
behaviour /t(— a;) ~ 1 — 7 — | In a;, where 7 ^ 0.577216 is 
Euler's gamma constant, therefore 



fr - 



1-7-1- In 2- In r. 



(A23) 



as shown in Fig. lAll Thus, it is a slowly varying function, 
and the potential, to good accuracy, is quadratic in t. 



APPENDIX B: BASIS OF SOLUTIONS FOR 
PROPAGATION EQUATION 

In this appendix we build a complete system of solutions of 
the equation 1341 1 by further separation of D and r variables: 
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Figure Al. Function /r I —-^ )' defined by JAISI (solid line) 
and approximation IA23I (dashed line). 



^'(AT)=^enUr) 



replacing the original problem with 



Ac 



(Bl) 



(B2) 



for the set of all possible A. For u ^ 0, which is of primary 
interest for this study, it is more convenient to divide both 
sides by —iuj/2c denoting A — —iujjj,/2c and get rid of k- 
proportional term by further transforming 



— i— k ■ T 



to obtain the following equation for Cm (''")• 



2c, 



c^2 



mC(t) = dlCi^) + i-V'{T,cu) + ^k' C(t 



(B3) 



(B4) 



Now, introducing cylindrical coordinates t = re^, and sep- 
arating variables again, while observing periodicity 



C^ii-rex) ^Cn^{r)e 



n £ Z, 



(B5) 



(B6) 



(B7) 



we arrive at an ordinary differential equation 

with 

v„{t) = i—V'{t,uj) + — k" - — 

It is clear that the term depending on k can be absorbed 
into fj., so for each k mode the resulting equation has the 

same form with i-n^ = fj, 7 k . We will see soon that for an 

oscillator-like potential additional constraints on ^n^ lead to 
a discrete set of /i like in the corresponding QM case. 

Neglecting fr term in V'{t) as explained in Section 3 
and Appendix EI one rewrites the potential as 



v„{t) = — ( — +iK(a;,r)r^ 



with 

iK(cj) 



3 exp[iuj/ujo] 



[uj/ujoY 



\-Kv m A. 



(B8) 



(B9) 



Note that iK(a;,r) = (ift:(— a;,T))*; this guarantees that the 
solution to the problem with real initial conditions remains 
real for D > 0. 

We thus arrive at the differential equation 



M.<Cn.(r) = i A (,^) _ (^^ + i,(^),2 j ^„,,(,); (BIO) 

we will drop k and u) indices for now. The differential part 
of the equation above is Hermitean with respect to the stan- 
dard weight in problems with cylindrical symmetry (drr). 
However, since the oscillator constant is not real, the total 
operator is not Hermitean leading to different modes decay- 
ing at different rates. Therefore, we also need to solve the 
equation conjugate to IBlOl l: 



r1 I -^ Id/ dCL 



CL(r) (Bll) 



solutions to which will form a conjugate set to the set of 
solutions of IBIOII : 



(Ci,CM) = o, 



/m- 



(B12) 



Our aim is to represent the initial conditions as a linear 
combination of the solutions to IBlOII in the most economic 
way; the coefficients of this linear combination are found 
using solutions to (IBIH . In fact, this 'non-Hermitecity' is 
easily eliminated by a simple change of the independent vari- 
able and the solutions to IBlOII and IBIH can be obtained 
from those of a single equation. Indeed, if Ur{z) equates 



rur{z) 



1 d 

z Az 



AUr{z) 



+ Z I Ur{z) 



(B13) 



then it is easy to check that CriT) = iij-(cr) and Cr(''") = 
Ur(c*r) are the solutions of l|Blfll IbTT^ with /j, = rc^ and 
a — r'lc*)^ = fi* , provided that 



c = IK. (B14) 

Here one rather technical but important observation 
concerning the choice of c is due. The problem IIB13II does 
have a compact set of solutions Ur{z) which can be used 
to represent any square-integrable on the real positive ray 
function g{t), and all r in this set are real. Therefore, given 
a function /(r) we can always represent it as a linear com- 
bination of Cp if /(t/c) is square-integrable on r € (0, -l-oo). 
In all cases of interest, the required convergence only takes 
place at Re z > half of the complex plane, therefore among 
the four roots of IIB14II we will choose ones with real pos- 
itive part (which of the two is relatively unimportant as 
long as it is used consistently, we will use the one closest to 
the real axis because such choice is generally less restrictive 
on /(t)). To distinguish the appropriate c we will use the 
square root notation c^ = \/Tk, square root being under- 
stood in the above, arithmetic, sense. Let us also note that 
for real r 1B12I I changes into 

(C^Cr')=0, r^r'. (B15) 

To solve the equation lB13t we first change the depen- 
dent variable 



u{z) — v±{z) 



±"exp 









(B16) 



and then a change of the independent variable v{z) = w{z'^) 
brings it to the confluent hypergeometric form: 
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twl{t) + (1 ± n - t)w'±{t) - (ii^ + w±(t) = 0, (B17) 

where prime denotes the derivative with respect to the func- 
tion w argument t. 

The two independent solutions of the last equation are 
given by the confluent hypergeometric functions of the first 
and second kinds: 

■ 1 ± n r , , \ 
+ -,i±n,i 



The solution to the original problem 13011291 is therefore 
given by 

M(D) = Y^ w{Lo,D)e^^'-"^'' {mLi^,Mo) M^^^i.. (B24) 



with the basis functions, their conjugates and corresponding 
eigenvalues being 



W±{t)^CMM(l 



~CuU(^ 



2 4' 



(B18) 



Mkn!,^ (r, x) 



(B25) 



, 1 ± n, t 



2 4 

Cm,Cu are complex numbers. 

The solutions to <B13ll are therefore 

2 



Tv{\n\ +iy. 



((i.)^/V)'"'Ll"l(Vi;:r^) 



X exp 



i / \ in 

U]^[(z) — z exp 



and 



Ujj(z) — z ^ exp 



Af(i|^ + ^,l±n,.^ 



,, / 1 ± n r , , 



-i I rix H fci" cos(x ^ 0) I ^ Vi^ 



(B19) 



(B20) 



A^L.(^,X) 



(B26) 



l\^- 



7r(|n!+0! 






Of course, only two of these four are linearly indepen- 
dent as can be readily established JAbramowitz fc StegunI 
|l972)- The most universal choice of the two is um(z) = 

M{^ + l,l + \n\,z'') and uu{z) = 



X exp 



-i [nx + -fcr cos(x - 0) j - V^ 



and 



^'"1 exp 



z'"! exp 



IC , 9 ILO 



2lu 



kV— ^iKH(l + lnl+2Z). 



(B27) 



U\ 



1 + |n|, 2^). Again, only one of Here k = ke^, r = re^ and w{uj, D) is defined by ^. 



2 '4' 

these two, Um is regular at t = while ujj has a pole of or- 
der \n\ here. Since we are not interested in singular solutions 
the solution of interest to us is 

2 



Ur{z) = Z' 



exp 



M 



1 + F , ^ 1 , I I 2 



(B21) 



Now we recall that Re y'ire(a;) = Re •^— iK*(tj) > and 
from the asymptotics of the conflu ent hypergeometric func- 
tion (JAbramowitz fc StegurJll972 ) it follows that ("r and Q. 
grow exponentially with r ^ oo unless the first argument of 
M is a real, non-positive, integer number — /, Z = 0, 1, 2, ... in 
which case M{—1, 1 -1- |n|, t) becomes a polynomial of degree 
I in t known as the associated Laguerre polynomial of order 
\n\ and degree I L\"'(t). 

These polynomials form a basis in £2(7?."'") and therefore 
can be used to represent any initial conditions M^^^{t) of 
interest. This solves the problem; the complete eigensystem 
of IB 1 Oi l within regular at zero £.2{'Tl'^) functions is 



APPENDIX C: APPROXIMATION FOR THE 
SPREAD FUNCTION 

As discussed in the section 4, function w{uj,D) defined 
by H32|l represents the Fourier transform of the 'spread func- 
tion' w{t — t' , D) with which one has to convolve the solu- 
tion of problem 13411 . This spread function is universal in the 
sense that it is the same for all angular modes and thus phys- 
ically describes variance in the potential time delay within 
the ensemble of deflector field realizations. 

The inverse Fourier transform of 'w{u>) simplifies in the 
case when parameters of the problem (such as v, R, m) 
are constant, so that w{uj,D) = exp[7rt/_D(i?^A/ij +m^f^)]. 
In addition, most of the variability we are concerned with 
is concentrated at u> range where ui/ujo <C 1 and 'w{u)) is 
approximately Gaussian, according to IA19llX20l : 



OW = Cirl"l 



exp 



IK 2 

T 



(B22) 



w{u!) ~ exp 



m^o;^ / 2 4 2,3, 

-TTuU -— H H — m m A 

8c^ V 3 



(CI) 



fii = -2Vi^(l + |n| +2/), 



Z = 0,l,2,. 



while the complex conjugates of the above function are the 
complete eigensystem of IBlll : Q — Q, ai = nl . 

It is convenient to normalize them to unity (C; , 0') = 
5\i, which implies 



The inverse transform of this function is again a Gaussian 



2l\ 



[\n\+l)\ 



(iK) 



(|n| + l)/2. 



w{t — t') — yl-Ka^ j exp 



with dispersion 
2,^N m^ ^^2 f 41n^ A 



2af 



(B23) 



(C2) 



(C3) 



here we used that (LJ" (t))^f" exp[— t] is regular between 
argi = argViK and argf — 0, the integral over the 
arc connecting these rays approaches zero as its radius 
goes to infinity as well as orthogonality relations from 
iGradshtevn fc RvzhikI Jl994l) . 



and given that the time it takes light to traverse gravita- 
tional radius of the deflector of mass m is just tm = m/2c, 
and the total number of deflectors is N{D) — nuR^D, the 
effective width of the spread function w is 



. ^N rrr I 4 In^ A 

CTt{D) = tmVN\l 1 + ^ U 



N 



(C4) 
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The result might seem somewhat surprising in a sense that 
this dispersion depends on the total number of deflectors 
regardless of their location with respect to the observer in 
the final plane. However, the analysis we develop in this 
work does not involve an observer up to the point when some 
interpretation of the result is needed. Indeed, the variance 
above is a quantity defined on the plane labelled D (or, more 
precisely, the ensemble average of such planes) and therefore 
one needs access to the entire plane in order to measure it. 
When the observer can sample only a fraction of the 
plane, the dispersion will be less; moreover, for uncorrelated 
deflectors it is generally proportional to the square root of 
the number of deflectors - but only those which influenced 
the rays which have been observed. Our conjecture is there- 
fore that in calculating ljC4fl one only needs the number of 
deflectors within the region sampled by the observer over 
the run of the data set used. Here we touch the question of 
how average quantities measured on a distribution sample 
relate to those of the distribution itself, which is a complex 
problem on its own studied in statistics. A natural way to 
quantify the uncertainty of such measurements is provided 
by considering their second moments. This is the subject of 
the companion paper. 



